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ABSTRACT 


Presence of cracks in hull plates of submar ine/shi p is 
undesirable as it reduces the strength and efficiency of the 
vessel. Further underwater shocks created by mines torpedo etc. 
might lead to propagation of these cracks. It is important to 
know the parameters which lead to propagation of these cracks, 
resulting in a catastrophic failure. Stress intensity factor 
(SIF) is the most widely used parameter in static problems. 

The problem of submarine hull modelling is complex as it 
falls under the framework of dynamic elasto-plastic analysis. As 
a first attempt one can do static, linear elastic analysis to get 
an insight into the complex nature of the problem. In the present 
work a part of submarine hull bounded by two stiffening rings is 
idealised a plate. It is assumed that cracks are either parallel 
or perpendicular to stiffening rings, and the loading is 
transverse and uniform. As the submarine hull plate is 
considerably thick, Hindi in's plate theory has been used in 
formulating the problem. The problem is solved using the Finite 
Element technique. Stress intensity factor in mode I has been 
evaluated using the J integral approach. 

Effect of plate length, distance between the edges 
perpendicular to stiffening rings, crack orientation and the 
boundary conditions on the edges perpendicular to stiffening rings 
on the normalised stress intensity factor have been studied. It 
has been found that SIF increases as the distance between the 
edges perpendicular to stiffening rings is increased It is also 
found to increase as the crack orientation changes from parallel 
to perpendicular, and the boundary conditions change from simply 
supported to free on edges perpendicular to stiffening rings. It 
is found that a perpendicular crack with free ends is the most 


dangerous . 



ABSTRACT 


Presence of cracks in hull plates of submarine/ship is 
undesirable as it reduces the strengtli and efficiency of the 
vessel. Further underwater shocks created by mines torpedo etc. 
might lead to propagation of these cracks. It is important to 
know the parameters which lead to propagation of these cracks, 
resulting in a catastrophic failure. Stress intensity factor 
(SIF) is the most wi4ely used parameter in static problems. 

The problem of submarine hull modelling is complex as it 
falls under the framework of dynamic elasto— plast i c analysis. As 
a first attempt one can do static, linear elastic analysis to get 
an insight into the complex nature of the problem. In the present 
work a part of submarine hull bounded by two stiffening rings is 
idealised a plate. It is assumed that cracks are either parallel 
or perpendicular to stiffening rings, and the loading is 
transverse and uniform. As the submarine hull plate is 
considerably thick, Mindlin’s plate theory has been used in 
formulating the problem. The problem is solved using the Finite 
Element technique. Stress intensity factor in mode I has been 
evaluated using the J integral approach. 

Effect of plate length, distance between the edges 
perpendicular to stiffening rings, crack orientation and the 
boundary conditions on the edges perpendicular to stiffening rings 
on the normalised stress intensity factor have been studied. It 
has been found that SIF increases as the distance between the 
edges perpendicular to stiffening rings is increased It is also 
found to increase as the crack orientation changes from parallel 
to perpendicular, and the boundary conditions change from simply 
supported to free on edges perpendicular to stiffening rings. It 
is found that a perpendicular crack with free ends is the most 


danoerous . 



LIST OF FIGURES 


Figure Title Page 

2.1 - 12 

2.2 - ia 

2.3 - 20 

2.4 - 24 

E.5 Reduced and Selective 29 

Integration Rules 

2.6 - 32 

3.1 Test Problem for Singular Elements 35 

3.2 Deflection of Simply Supported 36 

Cracked Plate (Reference result) 

3.3 Deflection of Simply Supported 37 

Cracked Plate (Result obtained) 

3.4 Test Problem (J Integral) 39 

3-5 Variation of normalised SIF 40 

for a Cracked Square Plate 
(Reference result) 

3.6 Variation of Normalised SIF 41 

for a Cracked Square Plate 

(Result obtained) 

3.7 Plate Geometry and Crack 42 

Orientation 

3.8 Effect of h/B Ratio on Normalised 45 

SIF (Simply Supported Ends, Crack 

Parallel to Fixed Supports) 

3.9 Effect of L/B Ratio on Normalised 46 

SIF (Simply Supported Ends, Crack 

Parallel to Fixed Supports) 

3.10 Effect of h/B Ratio on Normalised 47 

SIF (Simply Supported Ends, Crack 
Perpendicular to Fixed Supports) 

3.11 Effect of L/B Ratio on Normalised 48 

SIF (Simply Supported Ends, Crack 
Perpendicular to Fixed Supports) 

ii 



Page 


F i gur e 
3. 1Z 

3.13 

3. 14 

3. 15 

3.16 

3.17 

3.18 

3.19 

3.20 

3.21 

3.22 

3.23 

3.24 


Title 

Effect of Crack Orientation on 
Normalised SIF (Simply Supported 
Ends, L/B = 1.0, h/B = 0.025) 

Effect of Crack Orientation on 
Normalised SIF (Simply Supported 
Ends, L/B = 1.0, h/B = 0.20) 

Effect of Crack Orientation on 
Normalised SIF (Simply Supported 
Ends, L/B = 1.5, h/B = 0,025) 

Effect of Crack Orientation on 
Normalised SIF (Simply Supported 
Ends, L/B = 1.5, h/B = 0.20) 

Effect of Crack Orientation on 
Normalised SIF (Simply Supported 
Ends, L/B = 1.8, h/B = 0.025) 

Effect of h/B Ratio on Normalised 
SIF (Free Ends, Crack Parallel to 
Fixed Supports) 

Effect of L/B Ratio on Normalised 
SIF (Free Ends, Crack Parallel to 
Fixed Supports) 

Effect of h/B Ratio on Normalised 
SIF (Free Ends, Crack Perpendicular 
to Fixed Supports) 

Effect of L/B Ratio Normalised 
SIF (Free Ends, Crack Perpendicular 
to Fixed Supports) 

Effect of Crack Orientation on 
Normalised SIF (Free Ends, L/B = 1.0, 
h/B = 0.06) 

Effect of Crack Orientation on 
Normalised SIF (Free Ends, L/B = 1.0, 
h/B = 0.20) 

Effect of Crack Orientation on 
Normalised SIF (Free Ends, L/B = 1.5 
h/B = 0.06) 

Effect of Crack Orientation on 
Normalised SIF (Free Ends, L/B =1.5 
h/B = 0.20) 


49 


50 


51 


52 


53 


56 


57 


58 


59 


60 


61 


62 


63 


iii 



LIST OF SYMBOLS 


a 

B 

CB3 

CDDj 

CDD 


Cf 

1 

Cf 

1 

CF> 

h 

J , J 
X y 

CJ3 

CJD-' 


K^,K2,K3 

•*V 

*<1 

Oil 


L 

LUl 

P 

z 

r 

t 

U 

u 

V 

w 

U 


Semi crack length 
Half plate width 
Strain displacement matrix 
Elasticity matrix (bending) 

Elasticity matrix (shear) 

Modulus of elasticity 
Body force vector 

Surface force vector on specified boundary 

Global force vector 

Plate thickness 

Path independent integrals 

Jacobian matrix 

Inverse of CJJ 

Stress intensity factors in mode I, II and III respectively 

Normalised stress intensity factor 

Global stiffness matrix 

Length of plate 

Shape functions 

Uniform transverse load 

Radial distance from crack tip 

Plate thickness 

Strain energy density factor 

Displacement in x direction 

Displacement in y direction 

Displacement in z direction 

Strain energy density 


IV 



Strains 


Greek Svmb 
(c) 

& Rotation about y axis 

& Rotation about negative x axis 

y 

v> Poisson ratio 

K tff rK Natural coordinates 

Curvatures 
£<y> Stresses 

y Rotation about y axis 

Vy Rotation about x axis 


V 



CHAPTEAR I 


INTRODUCTION 

1.1 INTRODUCTI<»{ 

Submarine hull or pressure hull is basically fabricated by 
welding of plates which are mounted on circumferential stiffening 
rings. Presence of cracks in hull plates is undesirable as it 
reduces the strength and efficiency of the vessel. Cracks can be 
present as inherent defects in manufactured components or may 
originate during service due to damages. Cracks may also emanate 
from stress concentration points due to fatigue or may be present 
as defects in welded joints during fabrication. The capacity of a 
plate to withstand external pressure pulse/transverse load reduces 
due to the presence of cracks. Further high intensity underwater 
shocks created by mines torpedo etc. lead to propagation of these 
cracks which might lead to catastrophic failure. As the residual 
strength is a function of crack size, it is essential to analyse a 
cracked plate in bending so as to determine the extent to which 
strength is reduced, and how this is related to crack size. The 
analysis also helps in determining maximum size of crack that can 
be tolerated and, time required for the crack to grow to critical 
size. 

A hull plate mounted on two adjacent stiffening rings can be 
localised for analysis/study of crack propagation in submarine 
hull due to underwater explosion. As the pressure pulse created 
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by underwater explosion acts -for a short duration of time and as 
the submarine hull behaves like an elasto plastic material at high 
stress level, the problem falls under the frame-work of dynamic 
elasto-plasti c analysis. A dynamic problem with elasto plastic 
mataerial behaviour will correctly model the problem and will give 
accurate results for fracture parameters, but it is difficult to 
analyse such problem. A rough estimate of fracture parameters can 
be obtained by analysing the problem within the framework of 
static elastic analysis. Although the values obtained may be 
approximate, the trend of variation of fracture parameters with 
respect to geometric variables such as crack length, plate 
thickness, aspect ratio of plate etc., is similar to the one 
predicted by dynamic elasto-plastic analysis. Thus static elastic 
analysis gives an insight into the complex nature of the problem 
and forms a base for dynamic elasto plastic analysis. 

1.2 LITERATURE SURVEY 

Stresses in the region near the crack tip are large and their 
distribution is complex. Several investigators have discussed the 
nature of the local stresses around a sharp crack in plates 
subjected to transverse loading. Williams C1T found that elastic 
bending stresses near the tip of a semi-infinite crack vary as the 
inverse square-root of the radial distance from the crack front. 
The magnitude of the stresses was left undetermined in his 
analysis. It was Sih and Rice CZ,3T who later determined this 
magnitude by using the theory of complex functions. The results 
in C1— 33 were obtained from the classical fourth order Kirchoff 
theory of thin plates, which has two short comings. The first one 
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is that the edge conditions at the crack surface are satisfied 
only in an approximate manner in that the three physically natural 
boundary conditions of prescribing bending moment, twisting moment 
and transverse shear are replaced by two conditions. Owing to 
such a replacement, the stress distribution in the immediate 
neighbourhood of the crack edges will naturally be affected and 
will not be accurate. The second limitation of Kirchoff theory is 
that transverse shear effects are neglected. Near a crack tip in 
a plate, the gradients of the in— plane stress components are very 
large. This induces a large variation of the transverse 
displacement in the plane— of the plate. This means an extremely 
large transverse shear strain and the associated transverse shear 
stresses . 

The theories which include the effect of transverse shear 
deformation and which make it possible to satisfy all the boundary 
conditions along a stress free edge are due Mindlin C4D and 
Reissner C53. Knowles and Wang C63 applied the Reissner theory 
making use of the integral transform technique, to obtain an 
approximate solution for the crack tip stresses in a plate with 
vanishingly small thickness. It was found that the order of the 
singularity of the boundary stresses near a crack tip was the same 
as classical theory CH but the angular distributions were not 
dependent on Poisson' s ratio like in the classical theory. 
Hartrancraft and Sih C73 applied the Reissner theory using another 
method of solution to study the effect of plate thickness on the 
bending stresses in the vicinity of the crack— front. He showed 
that, unlike in classical plate theory, the stress intensity 
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factor (SIF) depends strongly on the plate thickness. For a 
typical problem, the SIF was found to vary very rapidly in the 
range of 0 < h/a < 0.S5 where h is the plate thickness and a is 
the half crack length. 

The above solutions based on the Reissner theory establish 
only the functional form of singular stresses and not the general 
crack tip solution of William's type LiJ. Such a general solution 
was first provided by Murthy et al C83. The solution includes 
complete class of solutions expressed in terms of polar 
coordinates at the crack tip. It was observed that standard 
numerical techniques like collocation, successive integration etc. 
fail to give satisfactory convergence for SIF in plates with large 
plate sizes (1/a > 1), and small thickness ratio (h/a < 0.5). It 
was shown that such limitations in numerical analysis arose mainly 
due to the mathematical nature of continuum solution which 
involves Bessel functions with large arguments. This is a severe 
limitation because the parameters 1/a and h/a, in practice, do 
fall in the ranges mentioned. 

Although continuum solutions do give accurate picture of the 
singular stresses, their applicability is confined mostly to 
infinite domains. For general applicability, the finite element 
method (FEM) appears to be a logical choice because of its 
elegance, flexibility and applicability to large classes of 
geometries, materials and loading conditions. The various finite 
element formulations, with special reference to plate, bending 
problem, have been reviewed by Uahaba C93. Following three 
approaches are usually followed in treating the singularity while 
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formulating special displacement methods. 

In the first method, general crack-tip solutions from a 
continuum solution are chosen as the displacement field for the 
crack tip element. For convenience (SIF) is chosen as a degree of 
freedom. However, very few such formulations are available for 
plate bending case. Most of them are based on classical thin 
plate theory and thus do not account for the effect of plate 
thickness. However, the special crack tip element developed by 
Viswanath et al C10D is based on the sixth order plate theory. As 
a result it satisfies all the three natural conditions on the 
crack surface and accounts for the transverse shear deformation. 
The element can be of any arbitrary shape. A special feature of 
the element formulation is that only line integrations are 
involved in setting up of element stiffness matrix. Numerical 
studies reveal that there exists, an optimum crack tip element 
size for the given plate thickness which gives the best possible 
value of the stress intensity factor. The other limitation of 
this method is that nodal displacements become rather inaccurate 
at an infinitesimal distance from the crack. This limitation 
however can be overcome by refining the Mesh near the crack tip. 

In the second method, only the term corresponding to singular 
component of stress is considered from the analytical solution and 
is superposed on regular polynomials for the displacements in 
element development. Yagwa and Nishioka C113 have used this 
method to develop an &-noded isoparametric element for plate 
bending. Although the effect of transverse shear deformation is 
included in the element, the normal displacement at the crack tip 
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is still based on the Kirchoff assumption. 

In the third alternative, singularity is achieved through the 
Jacobian determinin'^ isoparametric element vanishing at the 
crack tip. A number of special crack tip elements have been 
developed from the standard isoparametric element enabling the 
singularity of the strain field et the crack tip to be modelled. 
Yamada C12D has summarised the development of singularity elements 
and has elucidated the correlations between proposals which take 
seemingly different forms. He has also shown that significant 
improvement is Achieved in the estimate of stress intensity factor 
by the use of singular elements. Barosum C13,14l has developed 
the quarter point quadratic singular element achieving the inverse 
square root singularity characteristic. The singularity is 
achieved by placing the mid side nodes near the crack tip at 
quarter points. The element contains the necessary rigid body 
modes and constant strain modes. Further, the element satisfies 
the continuity requirements of the displacements as well as 
rotations. The element is applicable to thick as well as thin 
plates. Since the shear deformation is included and it is 
possible to satisfy eu -tbe three natural boundary conditions, 
accuracy of the element is very high when compared with other 
results C73 obtained using Reissner theory. This is in spite of 
the fact that the element does not have the exact trigonometric 
angular distribution of stresses around the crack. It is found 
that four or more elements are sufficient to achieve this 
variation. Barosum C13,143 determined that by collapsing one side 
of a quadrilateral to form a triangular element in two dimensional 
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problems and one face of a brick to form a prismatic element in 
three dimensional problemSr far better results are obtained. 
Barosum also showed that triangular elements possess the same 
singularity in all directions. 

In addition they can have either 1/-^ r singularity 
representing perfectly elastic behaviour, or 1/r singularity 
representing perfectly plastic behaviour. These are achieved by 
having either one-node or multiple independent nodes at crack tip 
respectively. The ease of using these elements, as they exist in 
many general purpose programs, as well as their high accuracy 
makes them very attractive for application in linear fracture 
mechanics problems of plate and shell bending. 

There are various parameters called as fracture parameters, 
which can predict the initiation and propagation of cracks. If 
these parameters for a given problem can be determined 
analytically, numerically or experimentally and if critical values 
of these parameters for the given material are known from 
experiments, premature failures can be avoided. Most of the above 
references consider stress intensity factor (K) as the fracture 
parameter. Linear elastic fracture mechanics is based on the 
premise that the intensity of stress or strain field surrounding a 
crack tip may be uniquely described in terms of K. Thus, K can be 
viewed as a parameter characterising the resistance of metals to 
both static and fatigue fracture. The physical significance of 
the stress intensity factor approach is that it can be used as a 
design parameter to determine the nature of fracture process or 
used to determine the residual life of a component subjected to 
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dynamic loading. Besides depending on the applied loading, 
geometry of the plate, crack size and location it also depends on 
the mode of fracture. There are three basic modes associated with 
Plate Bending problems. They are 
Mode 1 t bending mode 

Mode II ! twisting mode 

Mode III ! shear mode 

First two modes are characterised by moment intensity factor while 
shear mode is represented by shear force intensity factor. The 
bulk of fracture mechanics work has been devoted to mode I 

analysis. This is due partly to the simplicity of its application 
and partly to its severity on the component. Cracks which start 

in single mode may in fact become mixed mode later on during their 

useful life. Since the concept of K is based on linear elasticity 
it has limitation both on accuracy and interpretation when 
materials capable of plastic deformation are considered. This 
limitation can be overcome by the J integral approach which is 
discussed in next paragraph. 

If the crack tip stresses are not modelled accurately, 

evaluation of K becomes inaccurate. The fracture parameter which 

avoids evaluation of crack tip stresses is the vector of J 

integral. Sosa and Eischen C153 and Sosa and Hermann II16I1 have 

proposed path independent integrals J and J for cracked plates 

X y 

subjected to bending loads using Reissner and Mindlin plate 
theories. When the strain energies of bending and shear are 
separated all the three stress intensity factors can be calculated 


8 



from the integrals and J^. This approach makes analysis of 
mixed mode problems simpler. The most outstanding feature of the 
J integrals approach is its path independency. This allows 
calculation of linear, non linear elastic energy release rate and 
elasto-plastic work remote from the crack tip. This concept can 
be extended to non-linearly elastic or elasto plastic materials. 

In case where linear elastic fracture mechanics is applicable 
the fracture behaviour can be determined by using strain energy 
density criterion. Sih C17,183 has proposed that failure occurs 
when stationary value of the strain energy density factor U 
reaches some critical value which is determined experimentally. 
On physical grounds, material elements are capable of storing 
energy by experiencing volume change (dilatation) and shape change 
(distortion). Excessive dilatation may result in brittle fracture 
while excessive distortion may tend to yield the material. The 
proportion of energy absorbed in dilatation and distortion can be 
determined for various material properties at locations of the 
stationary values of strain energy density factor. The direction 
of crack initiation as predicted from the strain energy density 
theory depends on the intensity of 1/r energy field near the crack 
tip. 

1.3 Objective and Scope 

The domain considered is part of the hull bounded by 
stiffening rings. The limiting dimension in the other direction 
will be determined from parametric studies. As the thickness to 
radius ratio is small the domain can be idealised as a plate. 

The edges mounted on stiffening rings are considered as fixed 
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and the distance between these edges is assumed to be constant. 
The conditions on other two ends do not fall in standard 
categories. Therefore they are parametr i cal ly taken as free and 
simply supported. 

The formulation has been done within the frame work of linear 
elastic fracture mechanics. The analysis is limited to static 
loading only. As the thickness of the plate is not small, 
Hindi in*s plate theory has been used to formulate the problem as 
it accounts for transverse shear deformation effect. Weld joints 
for securing plates on the stiffening rings are either parallel or 
perpendicular to rings. As cracks originate from weld joints it 
is assumed that cracks are either parallel or perpendicular to 
rings. 

Finite element method (FEM) has been used to solve the 
problem. Collapsed quarter point elements have been used to model 
the crack tip while quadrilateral isoparametric elements are used 
in the rest of the domain. In the present study, J integral 
approach has been used to evaluate the stress intensity factor 
(K). Because of symmetry of geometry, boundary conditions and 
loading. Kg and are very small, and as such only has been 
evaluated. Parametric study has been done to investigate the 
effects of crack length, plate length, plate thickness and crack 
orientation on SIF. 

1.4 Plan of Thesis 

After the introduction. Chapter II presents finite element 
modelling and formulation of the problem. Chapter III includes 
test problems, results of parametric study and conclusion. 
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CHAPTER II 


Finite Element Formulation 

In this chapter finite element formulation of a cracked plate 
subjected to transverse load, using Mindlin’s plate theory is 
presented. It has been assumed while formulating the problem that 
the material is isotropic, stress strain relationship is linear 
and the deformation is small within the zone of loading. Singular 
elements have been used at the crack tip while eight noded 
isoparametric elements model the rest of domain. J integral 
approach is used to evaluate the stress intensity factor. 

2.1 Mindlin’s Plata Theory 

The main assumptions of mindlin plate theory ares 

(a) Stress component along the normal to the midplane of the 
plate is negligible. 

(b) Transverse displacement does not vary in the thickness 
direction. 

Cc) In plane displacements of the midplane are negligible (no 
membrane effect). 

(d) The normal to the midplane remains straight but not 
necessarily normal to it after deformation. 

In classical Kirchoff thin plate theory it is assumed that the 
normal to the midplane remains normal to it even after 
deformation. This means that the cross sections do not’ deform (no 
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For a plate bending problem it is more convenient to deal 
with stress resultants rather than the stresses. The stress 
resultants are defined as follows i 
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where t is the plate thickness. Substitution of (2-5) in (2.6) 

and (2.7) gives stress resultants in terms of d , 0 , w. 

X y 


n 


12 ( 


M = 

y 12( 


Et“ f X ^ y 1 

I anr " “ s/ J 
El r __y + V 1 

lay ax J 


( 2 . 8 ) 


H 

xy 


_§!_ r » + „ y 1 

K1-K.) *=< J 

3x = I; ] 




Et 


y 2( 1-*v 


T ( ®y Ip ) 


(2.9) 


According to eqn. (2.5), the shear stresses and are 

constant over the thickness. To take care of variation of r and 

T the expressions for Q and Q are multiplied by a factor 
zx y X 

called as shear correction factor. Eqn. (2.9) then becomes 
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(2.10) 


where is the shear correction factor. For subsequent finite 
element formulation it is convenient to express equations (2.8) 
and (2.10) in matrix form as follows s 
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or 


and 


or 




r M 1 

X 



1 

V 

0 

M 

V 


Et^ 




► 

1£(1-»^^) 

V 

1 

0 

J 




( 1-J^ ) 

« — ** 




■ 0 

0 


where 


6y. 

y 

'¥g~ 

X 

3^ 


se 

+ ^ 


cm = CD>. ixy 

o 


f Q 1 


■ 1 o' 


r © + ^ 1 

X 

« Et 

— s 



K 


1-H> ) 




Q 


0 1 


+ J" 

y J 




^ X Sk J 


€Q> = CD> C«> 
s s 


■Cli> = moment vector 

{D>j^ = Elasticity matrix (bending) 

■C;t> = curvature vector 

CQ> = shear force vector 

= Elasticity matrix (shear) 


tel 


Shear strain vector 


as : 


The equilibrium equations for a mindlin plate may be 


SQ SQ 

__js + .__y + p =0 

ax ay z 


> 

^x 


- Q 

^y X 


m = 0 

X 


xy + ^ 

"Six 


Q + m =0 

y y 


(E.11) 


(E.1E) 


written 


(E.13) 
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where P 

2 





distributed force in z direction 
distributed moment about y axis 


<= distributed mo ment about negative x— axis. 
These equations have to be supplemented by boundary 
Typical boundary conditions are discussed in the next 


conditions . 
paragraph . 


The boundary C of the domain (Fig. E.2) is divided in 3 parts 
C^, Cg and C^. On transverse displacement and rotations are 
specified as s 



on Cg shear force and two moment components are specified 

* * * 

^ ^ ^ X X y y 

where 


(2.14) 


(2.15) 


on Cg, the 


q = Q n + Q n 

^ XX y y 


m 

X 

•4 

ti 



M 

X 

M 1 
xy 


n 

X 

m 

K. y J 

1 

M 

xy 

n 

y 


n 

L y J 


specified boundary conditions are mixed s 

* * * 

^ vs: ^ f ifi — III 

' s s n n 


(2.16) 


(2.17) 


(2.18) 


where 


6 = e S + © S (twist) 

s XX y y 

(2.19) 


«n_ 


= ra 


n 


m n , 

y y 


(bending moment) 


In general the problem consisting of differential equations 
(2.13) and boundary conditions (2.14), (2.15) and (2.18) is 
difficult to solve analytically. The problem can be solved using 
the finite element method, for which the total potential energy 
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The Finite element method is the most widely used technique 
which finds an approximate solution by minimising the functional 
associated with given problem. For the present problem, the 
functional is the total potential energy given by equation (2.20). 
The method consists of dividing the domain into a finite number of 
elements assuming suitable interpolations for the unknown 
variables over each element and finding the unknown constants in 
interpolating functions by minimising the functional. 

2.2.1 Deriving finite element equations 

Finite element formulation based on Mindlin*s theory has one 
important advantage over the one based on classical thin plate 
theory. Mindlin plate elements require only continuity of the 
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lateral displacement w and the independent normal rotations 9 and 

X 

B . On the other hand elements based on classical Kirchoff thin 

y 

plate theory need C continuity of w. Thus Mindlin plate elements 

are simpler to formulate. In the present formulation 8— noded 

isoparametric elements are used CFig. E.3a). Compared to 9 noded 

Lagrangian elements, these elements are better as they avoid 

phenomenon of locking in the thin plates. 

The displacement and normal rotations at any point in a 

C 6 ) 

8-noded isoparametric element n are give by 



C2.E1) 


or 

8 

Cu> = E CN.3® fa. 
i = 1 ' ^ 

where N? are called shape functions, are known functions 
expressed in terms of the local coordinates and the 

displacement vector fa.>® contains unknown nodal values of w, & 

6 . For the sake of convenience rotations have been modified as 

y 

shown in Fig- E.3(b). Mathematically 



when a domain is divided into finite number of elements, the 
boundary enclosing the domain also gets discritised depending on 
the number of elements. Along a typical boundary element, 
displacement and normal rotations can be expressed as s 


a** . 
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0 0 
-N^ 0 

1 K 

0 -N 

I 


(E.E3) 


<u> = r CN.3*^ ia.y^ 
itl ^ ^ 

where represent one-dimensional shape functions. 
The curvature vector can be written as 


X + 2 



-Nf 

i,x 



(E.24) 


{»J = E CB3 £a > 
1=.1 

The shear strains can be written as 


Sw 

B 4- ^ 
X 


Sw 

& + ^ 

y 


(2.25) 


E CB3®. €a.>® 
i=1 ^ 


The rotation on boundary can be expressed as s 


& - & n + = E 

n XX y y 


L 0 - N-", - N^tiy J 


(2.26) 
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or 


^ J sx ^ ^ ^ ^ 


Substituting equations (2.22), (2.23), (2.24), (2.25) and 

( 2 . 2& ) in (2.20) we get the functional in terms of the unknown 


c h 

vectors fa.} and {a.} 

1 1 

. ”e 8 8 

” = ^ E E E CK^..3® fa.>® 

2 i = -j 1 bij j 

fit — 

e 8 8 


1“’"“ eT e p 

5 E E E J CK 3® fa - £ £ ia. 

A Slj J ^ 1 


ri #gk *111 m 

2 3 bp 3 3 hr 

-EE ff^ - E E if > 


Cf.>® 

1 


(2.27) 


b=1 i*1 


b=1 i=1 


where 


['^bij] 'I [®bi] K] [®bj] (bending stiffness matrix) 

^e ^ 

€ fit I © 

["sij] " -f Ki] K] [®.i] ■“ ttifr»..i Mlri.l 


{'.)■ --f [“.] 

^e 

{n}'^ =J H 


eT 


bT 


fbl dA 


{t> dS 


(body force vector) 


(2.28) 


=j K.]-": - 


(surface force vector on Cg) 


(surface force vector on C^) 


Here is the area of a typical plate element, Cj^ is the length 
of a typical boundary element, is the total number of plate 

elements, n.„ is the total number of boundary element on the part 

dH 

Cg and is the total number of boundary element on part Cg of 

the boundary. 
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Usintj a standard proced 


ure, the expression (E.E7) for the 
total potential energy is expressed in terms of the global 
quantities in the following manner 

n « ^ {A5 ' CKa {.Ay - fF> (2.29) 


Here CK3 is the global stiffness matrix which contains the sum of 
expanded versions of both and ^ over all the plate 

elements. Similarly the global force vector fF> contains the sum 
of expanded versions of all the three vectors ff.>®, Cf and 

bg 1 i “ 

■Cfi> over appropriate plate and boundary elements. The vector 

CA} is called global displacement vector which contains the nodal 

values of w, B and B at all the nodes. When the plate is in 
X y 

equilibrium, n takes the minimum value. Minimisation of n with 
respect to fAJ leads to the following algebraic equation 

CK3 fA> * {F> (2.30) 


At the crack lip stresses and strains are singular and they 
vary as inverse square root of the radial distance from the crack 
tip. If regular elements with biquadratic approximation for w, & 

X 

and are used to model the crack tip, then to get accurate 
stresses many such elements are required. There exist special 
crack tip elements which can model IZ-fr singularity of stresses 
and strains at the crack tip. Accurate results are obtained by 
using four such elements at the crack tip. By placing mid side 
nodes at quarter point (Fig- 2-4(a)), IZ-Tr singularity of the 
strains and stresses is obtained only along the sides 1-2-3 and 
1-8-7. If, additionally the side 1-8-7 is collapsed and node 8 is 
moved to the quarter point position (Fig. 2.4(b)), then this 
produces IZ-fr singularity in every direction. 
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Evftlu^liori of ift®t,r‘ix Anud: for^co vfiC'Lor' 


malr i ci*s 


The first step in calculation of the elemental stiffness 

•* r 1 ® , e 

^bijl i^dil evaluation of matrices 


j] and is evaluation of matrices [®bi] 

and involves differentiation of shape functions N® 
with respect to co-ordinates (x,y). This is done by using the 
chain rule : 


N® = N* „ ? + n! rt 

i»x i»Jr »x t fY) j X 


(2.31) 


N® « N® „ ? + N® n 

I ty 1 »y X,rt ,y 

The derivatives of (?»>)) coordinates with respect to (x,y) 
coordinates are evaluated from following transformation* 

X ««= x(? ,n) 

(2.32) 


y * y (? »n ) 


For isoparametric elements this transformation is given by 


r x 1 

6 

N? 

0 1 


e 

X . 

y J 

► * z 
1*1 

1 

0 

m 

1 


1 

e 


where (x?, y®) are the nodal values and N? are the same 
functions as used in the interpolation of w, and 6^ 

(2.21)), The Jacobian matrix of the transformation is given 


(2.33) 

shape 

(eqn. 

by 


CJ3 » 


■ Sk/SK 


Sy/SK ' 

SyZ&Tf 


(2.34) 
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p 8 

z 

1 1 

N®* , 
1 if 

e 

K . 

1 

6 

z 

i = 1 

N® „ 

1 »? 

8 

z 

L i^i 

N® 

1 ,n 

e 

K . 
i 

8 

E 

i_i 

N? y 
1 tW 

^ X ' 


f rt ^ 

and Ti 

can 

f X 

-y 

w X 

»y 


matr i x from as 
“ ^ . 


t) 


n 


= CJ3 

. I 

.y 

The integrals defining the elemental mat 


( 2 . 35 ) 


rice, 


and 


r T * “3 

1 sijj elemental vectors ^^i^ 

evaluated numerically by Gauss— Legendre numerical integration 

formulas. For that purpose, the variables of integration need to 


be transformed from (X: 

,y) 

to 

(f r» ) . 

When 

this 

is 

integrals become 









-e +1 

+ 1 


eT 






% 

0^ 

ti 

J 


CD3,. 

b 

hi] 

® det 

CJD dlf 

dy) 

-1 

-1 







— -.0 *^1 

hi] “-f 

+1 

r *1 

#T 



e 



J 

1 

m 

t 


CD3 

s 

hi] 

det 

CJT df 

dry 


■i -1 
+ 1 +1 


/f .1 * J J I^nJ Cb3 det CJ3 djf 


df) 


( 2 . 36 ) 


1 -1 
+ i f -,bt 


T ds 

Ct3' ^ dC 


w-j H 


The derivative ^ is evaluated by using the transformation between 
the global boundary coordinate s and the local boundary coordinate 
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c which is consistent with transformation (2.33). 

The requirement of only C® continuity of displacements and 

independent rotations in the Hindlin formulation allows the use of 

a wide range of interpolation schemes. However, when exact 

numerical integration scheme is used for thin plates with standard 

Mindlin finite elements, the stiffness matrix becomes very large 

leading to inaccurate results. This phenomenon, which is termed 

locking, is caused by the imposition oT the constraints y — y 

xz * yz 

= 0 in the shear strain energy terms for thin plates. The part of 
the stiffness matrix corresponding to transverse shear strain 
energy may be interpreted as penalty functions which force the 
shear strains to become zero as the plate thickness-to-span ratio 
is reduced. 

When a reduced order of numerical integration of the 
stiffness terms is used on some isoparametric Mindlin plate 
elements improved behaviour is often obtained. The use of reduced 
integration, though successful in some cases, produces the further 
problem of rank deficiency of the stiffness matrix. In selective 
integration approach reduced integration rule is used only for the 
stiffness matrix associated with the transverse shear strain 
energy. This is done to alleviate the overconstraining effects 
of this portion of stiffness matrix. Full integration is used on 
the remaining terms in an attempt to retain the required rank of 
the overall stiffness matrix. Fig. 2.5 shows the exact, reduced 
and selective integration rules used for typical serendipity and 
Lagrangian mindlin plate elements. 
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Glob.»l a*. ‘(PiTibly of elemental matrices .jj and j^K . .j 


to 


UK! and of element force vector Cf.>®, {f.> ^ and ff.J ^ into f F> , 

3. I i 

and application of essential boundary conditions on part and C„ 

1 3 

is done using standard procedures. 


2.2.31 Pout ProcftBssing and Strass Intanslty Factor 

The set of equations resulting after application of boundary 

conditions are solved using Frontal Solver. On solving these 

equations values of w, 0 and 0 are obtained. Further, 

X y 

curvatures, bending moment and shear forces are evaluated at gauss 
points . 


In plate bending problems when the distributed moments and 


m^ are absent certain Path independent integrals are used in the 
prediction of crack initiation. Uhen x-axis is taken along the 
crack, these integrals are defined by 


S0 
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00 
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f*^xy ^ ”y 0y 
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ir^ + M 

0K xy 

00 

y 
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( 2 . 37 ) 


^P. 


dO 


( 2 . 38 ) 




where 
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W (t w 
* z 


U strain energy density 

-■ distontinuity (or jump) across the crack surfaces, 
r - Path extending along the upper and lower crack 

surfaces. < 

r, ■ 

fi " domain enclosed by curve r 

n = outward drawn normal 


j can be broken into shear and bending part, such that 

X 

I I® j. 1* 

J *= J + J 

XXX 


(E.39) 


where 


d© ^ ^ ^ d©--j 

J® - J [w^ '' [^x JT * "xy 3^] "x ("xy THT * 3iin"yJ*®' 


y Wk 




(2.40) 


■"J [®X JT ' 

Ci 

39 ^ s 

K - -T "« " “y "y) 


iSCi 


The stress intensity factors are defined by the relations 



In the elastic 
are related to 



1 i» -fSr M < r ,0) 
r^O 

li® H^^(r,0) 

r-0 


lim Q (r,0) 

r-*-0 

case the stress intensity factors Kg and Kg 


j!’ K •'y “ 


iBi [" 1.2 . ,,2| 

ii? ri 

(2.42) 

6(1-H>)n ^2 

"3EK *^3 

(2.43) 

1 / |/ 

|\ M 

Eh 

(2.44) 
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on combining (2.42) and (2.43) 


J®+J^ 

X X 


lErr r 2 


Kg + ^ 


(1-K>) 


•< 1 ] 


(E.45> 


As the domain under consideration is subjected to Uniform 
transverse load and the geometry and boundary conditions are 
symmetric only is dominant and Kg, are negligible. 
Therefore is evaluated from by setting Kg and Kg to zero in 
equation {2. 45). In the present study crack is either parallel or 
perpendicular to the fixed supports. Since x-axis is along the 
crack in the definition of the direction of x axis in the two 
cases is different. The typical mesh and path of integration are 
shown in Fig. 2.6. 
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CHAPTER III 


Results & Discussions 

In this chapter the Hindlin finite element formulation 
developed in Chapter II has been applied to a part of submarine 
hull subjected to uniform transverse load* To check the software 
package results obtained by the same are compared with test 
problems at various stages. The variation of the stress intensity 
factor in mode I is studied with respect to the crack length, the 
plate length, the plate thickness, the crack orientation and the 
boundary conditions. 

3.1 Test Problems 

Two test problems discussed by Hinton and Owen C193 are 
solved to check the basic software for an uncracked plate in pure 
bending based on Mindlin's plate theory. The results obtained 
were in complete agreement with those presented in II193. 

3.1.1 Singular Elotnenis 

The problem of simply supported rectangular plate with a 
crack at the center is solved for the case of a uniform transverse 
load. Fig. 3.1 shows the domain and finite element model for a 
quarter of the plate. The total number of elements used in the 
model are 22 including four crack tip singular elements. The 
material properties of the plate are s 

E = 10920 V = 0.3 

mat 

Fig- 3.2 shows transverse displacements at different sections of 
the plate as presented in reference C133. Fig. 3.3 shows the 
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Fig- 3.2 


Deflection of a rectangular cracked plate 
(Reference result) 
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displacements obtained by the present model. It can be seen that 
the trends in the two figures match quite well. This shows that 
singular elements model behaviour of cracked plates quite well. 

3.1.2 J integral mrtd SZF 

To verify the validity of J integral approach to determine 
the stress intensity factor a test problem is solved and the 
results are compared with those of Sosa and Eischen C15T. The 


domain is a 

square 

plate , 

simply supported 

all 

around 

and 

subjected to a 

uniform 

pressure load, as shown in 

Fig . 

3.4. 

The 

finite element 

model contains 

66 elements in a 

quarter of 

the 


plate. The material properties used s 


E - 1000 w = 0.3 
mm^ 

Fig. 3.5 shows the result from Reference C15I1, on the effects of 
crack length and plate thickness on the normalised stress 
intensity factor. The normalised stress intensity factor is 
defined as 


K. 


P fi 
z 


where 

P « uniform transverse load 
z 

B » halt plate width 
a * semi crack length 

The results obtained by the present model are depicted in Fig. 
3.6. It is observed that the trends and the values are in close 
agreement with those shown in Fig. 3.5. This problem validates 
the complete formulation presented in Chapter II and the software 

package developed. 
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3. 2 Par Amftlr ic Si.udlac 

Submarine hull is fabricated by welding of plates on 
circumferential rings. One such plate mounted between two 
adjacent stiffening rings has been localised for analysis as 
shown in Fig. 3.7. The two ends welded on stiffening rings are 

idealised as fixed ends and the distance (2B) between them is 

assumed constant. As the boundary conditions at the other two 
edges do not fall in standard categories they are parametrically 
chosen as free and simply supported. The distance L between these 
two edges and the thickness h of the plate are varied to study 

their effects on the stress intensity factor. The crack is taken 

either parallel or perpendicular to the fixed supports (Fig. 3.7). 
The material properties used are : 

E = 1000 v = 0.3 


3. 2. 1 Simply Supported Case 
(A) Crack Parallel to fixed supports 

Normalised stress intensity factor is plotted as a function 
of a/B, for various h/B and L/B ratios in Fig. 3.8 and Fig. 3.9 
respectively. The normalised stress intensity factor is defined 
as before. In an infinite plate containing a crack of length 2a 
and bent by uniform bending movement applied at infinityr it is 

observed that the normalised SIF ia) increases with plate 

thickness h but decreases with the semi crack length *a 
(Hartrancraft and Sih Z71, Viswanath et al C103). In a finite 
plate also similar trends are observed by Viswanath et al C103 for 
the case of edge movements and by Sosa and Eischen C153 for the 
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case of uniform transverse load. This may be because the infinite 
plate solution for the bending moment distribution remains valid 
in the neighbourhood of a crack in finite plate as well. The 
trends exhibited by our results (Fig. 3.8) agree with those in 
above references. Note that although normalised SIF decreases with 
*a*, the dimensional SIF may not necessarily decrease with it as 
the expression for normalised SIF contains -fa in its denominator. 
Further it is observed that when the thickness to width ratio 
(h/B) is reduced beyond 0.02, there is hardly any change in the 
normalised SIF. In Kirchoff's thin plate theory, SIF does not 
depend on the plate thickness, so it can be concluded that the 
limiting h/B value upto which the thin plate theory is valid is 
approximately 0.02. 

As the length (L) is increased it is observed that normalised 
SIF increases (Fig- 3.9). The .stress intensity factor depends on 
the bending moment values in the vicinity of the crack. For a 
parallel crack these values are primarily affected by the 
reactions at fixed supports. However these reactions are 
influenced by the boundary conditions on other two edges as well. 
It is observed that beyond L/B = 2.0 normalised SIF is not 
significantly affected by these boundary conditions. 

(B) Crack perpendicular to fixed supports 

The trend of variation of normalised SIF with respect to h/B, 
L/B and a/B ratio is similar. However the values of normalised 
stress intensity factor are more in this case (Fig. 3.10 and Fig. 
3.11). This happens because when the crack is perpendicular, the 
bending moment values in the neighbourhood of the crack are not 
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influenced by the moment reactions at the fixed supports. Whereas 
the limiting h/B value, upto which the thin plate theory is valid 
is slightly more than case (A), on the other hand the L/B value 
beyond which effect of boundary conditions vanishes is the same. 

(C> Orientation effect 

It is seen from Fig. 3.12 that for a thin square plate (L/B = 
1.0, h/B = .025) a crack perpendicular to the fixed supports has 
the larger normalised SIF than crack parallel to supports. Thus a 
perpendicular crack is more dangerous than a parallel crack. As 
the thickness is increased (h/B == 0.20) the difference becomes 
smaller (Fig. 3.13). When the length is increased the difference 
becomes even smaller (Fig. 3.14). Fig. 3.15 shows that for a 
long and thick plate (L/B = 1.5, h/B = 0.20) the curves nearly 
coincide. Thus when both ’L* and ’h’ increase the normalised SIF 
becomes independent of the orientation. However for moderately 
thick plates the curves can cross each other when L/B is 
sufficiently large. This means that critical orientation now 
depends on a/B ratio (Fig. 3.16). It can be observed that his 
critical a/B ratio is around 0.45. 

3. 2. 2 Case with free ends 
(A) Crack parallel to fixed supports 

It can be seen from Fig. 3.17 and 3.18 that the trend of 
variation of normalised SIF with respect to the crack length, 
plate thickness and plate length are similar to the earlier case 
(i.e. when the edges are simply supported). When the other edges 
are free, the entire transverse load is supported by the fixed 
supports. This leads to substantially higher values of support 
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reactions which in turn raise the levels of bending moment in the 
neighbourhood of the crack to great extent. Because of this the 
normalised SIF values are significantly higher when the other 
edges are free. It is also observed that limiting values of L/B 
and h/B are also considerably higher (L/B = 2.2, h/B = .025). 

(B) Crack perpendicular to fixed supports 

When the other two edges are free, the case of perpendicular 
crack is very interesting. Here the trends of variation of the 
normalised SIF with respect to a/B, h/B and L/B are completely 
different than all the three previous cases. As far as variation 
with a/B is concerned, unlike in other cases it is quite sharp in 
the lower range of a/B. Further the change in the normalised SIF 
with respect to either the plate thickness or the plate length is 
very insignificant. It has been observed earlier that as the 
boundary conditions change from simply supported to free the value 
of normalised SIF increases. It is also seen that value of 
normalised SIF increases as the crack orientation changes from 
parallel to perpendicular. In this case, due to the cumulative 
effect of above two factors the value of normalised SIF is very 
high. 

(C) Orientation 

It has been observed in section 3.2.2(B} that for a 
perpendicular crack, the length and thickness have insignificant 
effect on the value of normalised SIF. So with increase in 'L* or 
•h* only the curve with parallel crack orientation shifts (Fig. 
3.21 — 3.24). As the variation with a/B is different, the curves 
cannot coincide completely. It can be seen that for lower values 
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of a/B the orientation effect is observed while for 
ratios the effect vanishes. 


higher a/B 


3.3 Coruclusiortc 


In the present work a part of submarine hull subjected to a 
transverse load has been idealised as a plate. It has been 
assumed that the two ends of the plate mounted on circumferential 
stiffening rings have fixed boundary conditions and the distance 
between them remains constant. Effect of plate thickness, plate 
length, crack orientation and boundary conditions on normalised 
SIF has been studied. The following conclusions can be drawn from 
this work. 


(1) It is well known that as the plate thickness (h) is 
decreased beyond a certain value, the thin plate theory 
becomes valid and as a result the normalised SIF becomes 
independent of it. In our case the limiting value of h/B is 
between 0.02 to .025, the lower value being applicable to the 
case of parallel crack and simply supported boundary 
conditions. 


(2> It is observed that normalised SIF increases with L/B ratio 
except for the case of perpendicular crack with free 
supports. It is further observed that for a simply supported 
case the limiting value of L/B ratio after which the effect 
of boundary conditions vanishes is 2.0. For the free case 
with parallel crack the limiting value is 2.2. 

(3> It is observed that for a simply supported case for L/B = 1.5 
and h/B == 0.20, SIF becomes independent of orientation, while 
for a free case such a behaviour is not observed. 
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(4) It is that value of normalised ^TP 

'-'rmaiisea SIF increases as the 

boundary conditions change fmm - . 

« trom simply supported to free, and 

as the orientation changes from n^r^iioi ♦ , 

on» parallel to perpendicular. 

It can be said that a peroena, --..t « • 

perpendicular crack with free ends is 

the most dangerous. 


3.4 Suggestloris; for future work 

The finite element »„d,l „„ improved by using e 
degenerate solid element «th five degrees per node, so that the 
membrane effect can be modelled accurately. The analysis can be 
extended to include elasto-plastic dynamic behaviour. The present 
idealisation of submarine hull as a plate can be modified to a 
cylindrical ring (shell) analysis to avoid uncertainty of 
boundary conditions- 
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